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We present a new test of gravitational interactions at the r ~ (0.2 — 20) Mpc scale, around the 
virial radius of dark matter halos measured through cluster- galaxy lensing of maxBCG clusters from 
the Sloan Digital Sky Survey (SDSS). We employ predictions from self-consistent simulations of f{R) 
gravity to find an upper bound on the background field amplitude of |/flo| < 3.5 x 10~^ at the ID- 
marginalized 95% confidence level. As a model-independent assessment of the constraining power of 
cluster profiles measured through weak gravitational lensing, we also constrain the amplitude Fq of a 
phenomenological modification based on the profile enhancement induced by f{R) gravity when not 
including effects from the increased cluster abundance in f{R)- In both scenarios, dark- matter-only 
simulations of the concordance model corresponding to |/ho| ~ and Fo — are consistent with 
the lensing measurements, i.e., at the 68% and 95% confidence level, respectively. 



I. INTRODUCTION 

Modifications of gravity can serve as an alternative ex- 
planation to the dark energy paradigm for the late-time 
accelerated expansion of our Universe. Such modifica- 
tions have extensively been tested on solar-system scales 
(see, e.g., [l|) and to a lesser degree at large cosmolog- 
ical scales using specific alternative theories of gravity 
(see, e.g., 0-0), as well as generic modifications to gen- 
eral relati vity (G R) while adopting a ACDM background 
(see, e.g., [l0l - [2l| ') or simultaneously allowing a dynamic 
effective dark energy equation of state [13, [l^]. However, 
gravity may also be tested by the structure observed at 
intermediate scales [H, H^. In this regime, nonlinear 
gravitational interactions gain in importance and need 
to be modeled correctly to obtain reliable predictions for 
both GR and its competitors, which in turn can be com- 
pared with observations to infer constraints on modified 
gravity theories. 

To study nonlinear effects in structure formation, we 
need to specialize to a particular gravitational modifica- 
tion. In our case, this is f{R) gravity. Within this model, 
the Einstein-Hilbert action is supplemented with a free 
function /(i?) of the Ricci scalar R. It has been shown 
that such models can reproduce the late-time acceler- 
ated expansion of the Universe without invoking dark 
energy |26l - [28j . However, they also produce a stronger 
gravitational coupling and enhance the growth of struc- 
ture. f{R) gravity is formally equivalent to a scalar- 
tensor theory where the additional degree of freedom is 
described by the scalaron field fa = d//di? [i^,!!^]- We 
parametrize our models by the background value of the 



scalaron field today, |/ijo|- The fn field is massive, and 
below its Compton wavelength, it enhances gravitational 
forces by a factor of 4/3. Due to the density dependence 
of the scalaron's mass, viable f{R) gravity models ex- 
perience a mechanism dubbed the chameleon effect 
|33| . which returns gravitational forces to the standard 
relations in high-density regions, making them compat- 
ible with solar-system tests [11] at r < 20 AU. The 
transition required to interpolate between the low curva- 
ture of the large-scale structure and the high curvature 
of the galactic halo sets the currently strongest bound on 
the background field, |/iio| < |*| (10^^ - IQ-^) [H, 
i.e., the typical depth of cosmological potential wells. A 
bound of the same order is obtained from galaxies serv- 
ing as strong gravitational lenses at r '--^ (1 — 10) kpc. 
Independently, strong constraints can also be inferred 
from the large-scale structure (r > 10 Mpc). The en- 
hanced growth of structure observed in f{R) gravity 
models manifests itself on the largest scales of the cosmic 
microwave background (CMB) temperature anisotropy 
power spectrum [ssj . where compatibility with CMB 
data places an upper bound on |//fo| of order unity 
Cross correlations of the CMB temperature field with 
foreground galaxies tighten this constraint by an or- 
der of magnitude [B-d, H^- However, the currently 
strongest constraints on f{R) gravity models from large- 
scale structures are inferred from the analysis of the 
abundance of clusters, yielding an improvement over the 
CMB constraints of nearly four orders of magnitude [3,[§1- 

In this paper, we present a new test of gravity at the 
r ^ (0.2 — 20) Mpc scale, i.e., around the virial radius of 
dark matter halos measured through the excess surface 
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mass density from cluster-galaxy lensing. A'^-body sim- 
ulations of modified gravity scenarios have shown that 
halo-density profiles exhibit a characteristic enhancement 
at a few virial radii when compared to halo profiles in GR 
simulations with the same expansion history ^36i . ISTj . In 
models which attempt to explain the accelerated expan- 
sion of the Universe without dark energy, the modifica- 
tions to the gravitational force generally increase towards 
late times, leading to a pileup of matter in the infall re- 
gions of massive halos. In contrast, the inner profiles of 
halos are less affected since they formed earlier, when the 
force modifications were weak or absent. 

Here, we use this effect to constrain the field ampli- 
tude Ifm] of the Hu-Sawicki (s^l /(R) gravity model 
with measurements of weak lensing around the maxBCG 
galaxy cluster sample [s^ from the Sloan Digital Sky 
Survey (SDSS) Through matching clusters by abun- 
dance, wc consistently take into account the modified 
gravity effects on the halo-mass function as well as the 
profiles. In addition, we consider a phenomenological ap- 
proach modeled on the f{R) effects on the halo profiles 
at fixed mass. While this approach is not entirely con- 
sistent (since it does not include the effects on the halo- 
mass function), the constraints are largely independent 
of halo number counts, and, moreover, are given directly 
in terms of the observable, rather than a model param- 
eter. They can thus be used to assess the constraining 
power of halo profiles measured through weak lensing on 
a wider range of modified gravity models, including, for 
example, models where gravity is weakened and profiles 
are consequently suppressed with respect to GR. In both 
cases, we perform a Markov chain Monte Carlo (MCMC) 
likelihood analysis on the underlying parameter spaces. 

The outline of the paper is as follows. In fjll we review 
the f{R) gravity model and weak gravitational lensing. 
We then describe the iV-body simulations employed to 
derive the dark matter halo properties ( pillj) . and the pro- 
cedure used to predict weak-lensing obscrvablcs in f{R) 
and ACDM cosmologies ( ^IV[) . ifVl then introduces the 
observational data as well as external priors used in this 
study. The constraints on the alternative gravity models 
are presented in Wli along with a discussion of system- 
atic effects that may contaminate the data or complicate 
its interpretation. We conclude in Will The appendices 
give further details about the halo model and interpola- 
tion used in mVl 



II. MODIFIED GRAVITY & GRAVITATIONAL 
LENSING 

When gravitational interactions are modified, the 
growth of structure and thus the distribution of mass, 
as well as the relation between light deflection and mass 
distribution change j40l - l43j . Effects of modified gravity 
on halo properties were studied in the case of f{R) grav- 
ity in, e.g., [3^ (cf. psj ) and the DGP brancworld 
scenario in, e.g.. (37l [46j (cf. (47|). 



We concentrate on Hu-Sawicki [SJI f{R) gravity and 
rely on the nonlinear behavior measured in A^-body sim- 
ulations of this model [Uliilii] (cf. [13, [EH). We shall 
first review the details of the Hu-Sawicki model and how 
to relate lensing observables to the underlying matter 
distribution. We then briefly review how stacked weak- 
lensing observables measure the mass distribution around 
halos. 



A. f{R) gravity 

In f{R) gravity, the Einstein- Hilbert action is supple- 
mented by a free function of the Ricci scalar R, 



S = 



J d*x^[R + fiR)]+ J d^x^C^. (f 



Here, Cm is the matter Lagrangian and wc have set c = 
1. Variation with respect to the metric g^i, yields the 
modified Einstein equations for metric f{R) gravity. 



Guv + fRRau 



f 



(2) 

where the connection is of Levi-Civita type and fn = 
df/dR is the additional scalar degree of freedom of the 
model, characterizing the force modifications. 

We specialize our considerations to the functional 
form [34| 



fiR) 
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C2 (R/m^y 



1 



(3) 



where = ^-KGpm/i- The free parameters of 

the model ci, C2, and n can be chosen to reproduce 
the ACDM expansion history and satisfy solar-system 
tests [131 through the chameleon mechanism jsil - lssj . In 
the high-curvature regime, cJ^^^R 3> rv? ^ Eq. ([3]) simpli- 
fies to 



ci 2 fmRo 
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f{R) = --m' 

C2 



(4) 



where Rq denotes the background curvature today, Rq 
R\z=o , and fno = fniRo)- We further infer 
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from requiring equivalence with ACDM when |/_ro| 
and restrict to models with n = 1. Varying n changes the 
evolution of the Conipton wavelength of the fn field with 
redshift. Generally, constraints on fpo become weaker 
(stronger) for n > I (71 < f ) (see [54I for a study of the 
mass function of halos in f{R) with varying n). In the 
following, we will further assume that |/i?,o| ^ 1, and 
drop terms that are higher order in fn. 

In the quasistatic limit, the trace and time-time com- 
ponent of the modified Einstein equations yield the fa 
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field equation and Poisson equation for tlie Newtonian 
potential ^ = i5(7oo/(2<?oo) in the longitudinal gauge. 
Specifically, 



a 



^2 



a 5p,, 
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■SRih 



(6) 
(7) 



Here, coordinates are comoving, S/b. ~ IbXR) — fB.{R), 
SR. = R — R, Spin = Pm — Pin- In contrast, the potential 
^_ = — $)/2, where $ = 5gii/{2gii), governing the 
propagation of light and hence lensing is only affected at 
order fa, which is of order 10~^ or less for the models we 
consider here. Hence, the modifications to cluster-galaxy 
lensing studied in this paper are caused by modifications 
in the distribution of matter, which arise from the en- 
hanced gravitational forces. 

If the background field \fm\ is large compared to typ- 
ical gravitational potentials 10~^), we may linearize 
the field equations via the approximation 



5R'. 



dR 



SfE 



(8) 



R=R 



where Ac = l/nif^ is the Compton wavelength of the 
field at the background. In Fourier space, the solution to 
Eqs. © and ([7]) within the linearized approximation is 



fc2^r(k) 




An — 



a^(5p,„(k), 
(9) 

where k = |k|. For scales k ^ 27rAQ^a, this leads to an 
enhancement of gravitational forces by a factor of 4/3. 
Computations using Eq. ^ are referred to as the no- 
chameleon or linearized f{R) case (49j . 

If the background field becomes small compared to the 
depth of the gravitational potential of the object con- 
sidered (|/i?o| ^ 10^^, small-field limit), the chameleon 
mechanism becomes active, suppressing non-Newtonian 
forces. More precisely, S/r fa —Jb and from Eq. ([6|), 
SR ~ SnGSpm, which restores the standard Poisson 
equation in Eq. ([7]). Given that the constraints on |/ijo| 
expected from our lensing data are well within the large- 
field regime (|//jo| ^ 10~^), we can apply the approxi- 
mation Eq. ^ in the simulations. 



B. Weak Gravitational Lensing 

Weak gravitational lensing serves as a powerful probe 
of the total matter distribution (baryonic -I- dark mat- 
ter) within our Universe. Here, we focus on stacked 
cluster-galaxy lensing, which measures the average defor- 
mation of background galaxy images around foreground 
maxBCG galaxy clusters. By averaging over many lenses, 
the contribution of unassociatcd large-scale structure is 
suppressed. 



We use the tangential shear 74 , measured using the el- 
lipticities of galaxy shapes, as a function of the comoving 
transverse separation from the lens r^ j w 9 (1 + z\) Di. 
Here, Di is the angular diameter distance to the lens and 
zi is the lens redshift. After stacking many clusters, the 
mass distribution becomes symmetric around the line of 
sight. Then, the shear is related to the excess surface 
mass density around the dark matter halos hosting the 
clusters, AT,{r±), through [53| 



lt{r±) 



AS](ri) 



-'crit 



(10) 



The excess surface mass density is related to the pro- 
jected surface density S(r^) through 



AS(ri)=E(ri)-S(ri) 
2 



E(rl)rldrl. 



(11) 



Here and throughout the paper, r denotes a three- 
dimensional separation, while r± refers to a projected 
two-dimensional separation. The comoving critical sur- 
face mass density is given by 



AnGDMl + zi)^ 



(12) 



where Dg and Dis denote the angular diameter distance 
to the source and between the lens and the source, respec- 
tively. Note that both Scrit and the conversion between 9 
and r± are dependent on the specific cosmological model 
(see WTCl . 

Assuming perfect centering of the lenses, the projected 
surface mass density is related to the halo profiles by 



S(r^) = 



4ttG 



Xl 



1 + Chr, 



.91 (x/ + y) 



r 



dy, (13) 



where H indicates the Hubble parameter, xi denotes the 
comoving distance to the lens, and y denotes the dis- 
tance from the lens along the line of sight. 6im(?') is 
the halo-matter correlation function which quantifies the 
total mass distribution around halo centers (see ^III Cp . 
The lensing window gi{x) depends on the source redshift 
distribution ps(x) as 



91 (X) = 2 



a[X)Ds{x) 



(14) 



assuming that Ps{x) is normalized to integrate to unity. 
The halo-matter correlation function decays strongly 
with increasing separation, so for the transverse scales 
considered in this study, the lensing strength gi{x) is ef- 
fectively a constant. 
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ibox [h ^ Mpc] Number of runs 
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128 


30 






64 


28 
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128 


6 






64 


6 




= 10"^ 


128 


30 






64 


28 


/ho 


= 10"^ 


128 


30 






64 


28 



TABLE I: Summary of f{R) simulation runs. All simulations 
use the linearized fn field equation, Eq. ((9|. The cosmological 
parameters of the simulations are given in ijlll Al 



III. SIMULATIONS 
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TABLE II: Parameter values for the zhorizon simulations: 
total and baryonic matter density parameters Qm and flh, re- 
spectively, the dimensionless Hubble parameter h, the power 
spectrum normalization erg — ag'^^'^, and the primordial 
spectral index Us- The first row indicates the fiducial cosmo- 
logical parameters inspired by the three-year WMAP best-fit 
values [5^, [s^] . 



In our study, we consider gravitational Icnsing mea- 
surements on scales of 0.5h~^ Mpc < r± < 25h~^ Mpc. 
These scales arc affected by nonlinear clustering such that 
numerical simulations are required to obtain reliable pre- 
dictions for the mass distribution. We utilize f{R) grav- 
ity simulations to obtain the deviations induced by the 
modified forces in the halo profiles with respect to the 
ACDM predictions, i.e., \fRo\ — 0, from the same initial 
conditions and simulation setup. We then employ the 
Ziirich Horizon (zhorizon) simulations which pro- 
vide ACDM predictions of better resolution and larger 
volume, and scale these results with the deviations from 
the f{R) gravity simulations. Note that we use simula- 
tions where the matter density field consists exclusively of 
dark matter particles, hereafter dark-matter-only (DM0) 
simulations. 



A. f{R) gravity simulations 

Since our constraints lie in a regime where the 
chameleon mechanism is not active and we require suffi- 
cient halo statistics, we employ no-chameleon f{R) grav- 
ity simulationSjWhich solve the linearized //j field equa- 
tion, Eq. ^ [3^, [4^, Simulations are conducted 
for l/flol = 10-^ 10-3, 10-^ and n = 1. Note 
that |//fo| = corresponds to ACDM. Other cosmologi- 
cal parameters are fixed to values following the WMAP 
3-year results, f^A = 0.76, f^b = 0.04181, h = 0.73, 
Us = 0.958, and the initial power in curvature fluctua- 
tions As = (4.89 X 10-5)2 at fc = 0.05 Mpc-\ corre- 
sponding to fJs — 0.82 at z = 0. The simulations are 
carried out on 512"^ grid cells with a total of iVp ~ 256"^ 
particles. Due to the limited volume and resolution of 
the f{R) simulations, we combine results from two dif- 
ferent box sizes, Lbox = 64/i-^ Mpc, 128h^^ Mpc. Only 
the smaller boxes contribute for r < 0.75h~^ Mpc, cor- 
responding to 3 grid cells for Lbox = 128/i"^ Mpc. The 
box sizes and number of runs for each value of \fRo\ arc 
summarized in Table [D 



Halos within the simulation and their associated 
masses are identified via a spherical overdensity (SO) al- 
gorithm (cf. HH). The particles are placed on the grid by 
a cloud-in-cell interpolation and counted within a grow- 
ing sphere around the center of mass until the required 
overdensity is reached. The mass of the halo is then de- 
fined by the sum of the particle masses contained in the 
sphere. This process is started at the highest overdensity 
grid point and hierarchically continued to lower overden- 
sity grid points until all halos are identified. The halos 
employed for this analysis (logj^Q Af > lO^^/i^^ -^o) gen- 
erally contain more than 10"^ particles. 



B. Concordance model simulations 

The ZHORIZON simulations comprise 30 -I- 24 pure dis- 
sipationless dark matter TV-body simulations of differ- 
ent ACDM cosmologies (see Table HI)) . designed for high- 
precision studies of cosmological structures on scales of 
up to a few lOO/i-i Mpc [54,, M]- 

The matter density field is sampled by A^p = 750'^ dark 
matter particles of mass M^m = 5.55 x 10^^h~^ Mq, in 
the fiducial case, with a box size of 1.5/1-1 Gpc. For the 
nonlinear gravitational evolution of the equal-mass par- 
ticles, the publicly available GADGET-2 code [s^ is used. 
In order to avoid two-particle collisions, a force softening 
length of 60h~^ kpc is employed. The transfer function at 
rcdshift z = is generated using CMBFAST (gO] and then 
rescaled to the initial redshift Zj = 50, where a realization 
of the potential on the grid is calculated. The particles 
are placed on a Cartesian grid of spacing Ax = 2h~^ Mpc 
and then displaced according to second-order Lagrangian 
perturbation theory using the 2lpt code [6l|, H^j • 

For each cosmology, we use four boxes from the 
ZHORIZON simulations, yielding an effective volume of 
I3.5h~^ Gpc'^. For all snapshots of each simulation, grav- 
itationally bound structures are identified by a Friends- 
of-Friends (FoF) algorithm [6^ with linking length of 0.2 
times the mean interparticle spacing |116l | . The halo cen- 
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ter is associated with the minimum of the potential of 
the particle distribution. Halos with fewer than 20 par- 
ticles are rejected, resulting in a halo-mass resolution of 
M > 1.2 X lO"'^'^ Mq, corresponding to a halo-number 
density n = 3.7 x IQ-^h'^ Mpc^ 



C. Cluster density profiles and sample selection 

Cluster-galaxy Icnsing measures a projection of the 
halo-matter cross correlation ^hm('')- We measure ^hm 
by averaging the spherically averaged density distribu- 
tion around halos in the ACDM and f{R) simulations: 



Pm 



- 1 



(15) 



In observations, clusters are selected according to their 
optical richness. The true mass can, however, deviate 
from the mass inferred from the mass-richness relation 
[gJ]- Thus, it is important to take into account the scat- 
ter in the mass-richness relation. In the simulations, we 
model the scatter by a log-normal distribution, assigning 
a new mass to each halo in the simulations by 



M 



exp 



ln(Mo)+A/'(0,CT) 



r' 



(16) 



where Mq is the true mass and Af is the normal distri- 
bution with zero mean and variance a^. The scatter a 
is left as a free parameter in the likelihood analysis. We 
apply this scatter, with a = 0,0.4,0.6,0.8, to the halo 
masses in the simulations, then we mass order the halos 
according to the simulated mass with scatter, and finally 
select the A^h most massive ones until the required clus- 
ter abundance n is achieved. Here, n is the estimated 
true average number density of the maxBCG sample (see 
Specifically, we require A^h ~ nVtot at z = 0.23, the 
mean rcdshift of the lens sample (see fV|) . where Vtot is 
the simulation volume. In the following, we will denote 
the corresponding mass profile as Chm^'^- The same pro- 
cedure is applied to the f{R) simulations at z = 0.22, 
but only for values of cr = 0, 0.6 for the scatter. 

As discussed in S|Vl a cylindrical cut is applied to the 
observational data in order to remove satellite galaxies. 
In the data this procedure removes both true clusters 
and satellite galaxies. To account for the removal of 
clusters, we mimic this approach in the ZHORIZON sim- 
ulation analysis, following the same algorithm and us- 
ing Sx = ±100/i~^ Mpc for the length of the cylin- 
der. This reduces the number density by 20% from 
n = 1.8 X IQ-^h-^ Mpc^ to n 1.45 x IQ-^h'^ Mpc^ for 
zero scatter and to n = 1.43 x 10~^h~^ Mpc'^ for scat- 
ter a = 0.4. Since the simulations contain only true halo 
centers, we conclude that 2/3 of the 30% of the maxBCG 
sample removed from the data were true clusters and 1/3 
were contaminating satellite galaxies. After applying the 
cylindrical cuts, the abundances of halos in the simula- 
tion and maxBCGs in the data sample agree very well. 



IV. FROM SIMULATIONS TO OBSERVABLES 

In this section, we describe how we obtain cluster- 
galaxy lensing predictions for /(i?) gravity from the sim- 
ulations described in the previous section. We also in- 
troduce our phenomenological approach modeled on the 
effects on the halo profile from f{R) modifications when 
averaging halos with the same lower mass threshold as in 
the concordance model. The intention of this approach, 
being largely unaffected by differences in halo number 
counts, is to yield a model- independent assessment of the 
constraining power of cluster density profiles measured 
through weak gravitational lensing. For this purpose, it 
is essential to not only study the f{R) modification on 
the abundance-matched halo profile but also its counter- 
part in a fixed mass range scenario as described in detail 
in gVH 



A. f{R) gravity halo profile predictions 

Since the f{R) simulations arc of worse resolution and 
smaller volume compared to the ACDM simulations, we 
parametrize the relative effect on the halo-matter cross 
correlation ^hm, rather than ^hm itself. That is, we mea- 
sure 



281111(7-, l/flol) 



Chin, Sim (?-, l/flol) 
Chm,sim(r, l/flol = 0) 



1 



(17) 



from the simulation outputs at z = 0.22 with |/ijo| > 
and l/flol = 0- We apply the scatter in mass as is done in 
the ACDM simulations (but only for a = 0, 0.6). In or- 
der to compare the f{R) gravity profiles to their ACDM 
counterparts, we consider two cases: a fixed common 
lower mass limit Mq, derived from the ACDM concor- 
dance cosmology (threshold- matched case, TM); and a 
lower mass limit for f{R) adjusted to match the abun- 
dance of tracers n (abundance-matched case. AM). Since 
the mass function of halos is enhanced in f{R) gravity, 
the f{R) mass threshold is higher in the second case. The 
AM case is a consistent approach for comparing f{R) 
gravity to ACDM; on the other hand, in the TM ap- 
proach, we purely rely on the modified gravity effects 
on halo profiles, without explicitly using the information 
from the mass function that has been used to place con- 
straints on f{R) in Q. 

The effects of a modification of gravity are significantly 
less severe in the TM case as compared to the abundance- 
matched case, i.e., when taking into account that mas- 
sive halos are more abundant in f{R). This effect is illus- 
trated in Fig.[Tl which shows Qsim(^) normalized to unity 
at the peak, i.e., g{r) [see Eq. p^ ]. and Fig. [2J which 
shows the peak amplitude as function of | /^o | • The pro- 
file enhancements peak at a few virial radii, correspond- 
ing to the infall region onto massive clusters. This ef- 
fect has also been found in simulations of other modified 
gravity models |37l |. and is a generic result of modified 



gravitational forces increasing towards late times (which 
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FIG. 1: The shape g(r) of the relative enhancement of ^hm('') in jiK) gravity simulations, for |/i?,o| = 10~^ in the abundance- 
(left panel) and threshold-matched case with scatter cr = (middle panel) and a = 0.6 (right panel), respectively. The middle 
panel shows the best-fit Gaussian function, Eq. (|20|l . to the simulation output. 



1.2 


AM, sim — ^ 






AM, model 




1.0 


TM, sim — ^ 




0.8 


TM, model 




0.6 






0.4 




--'"f" - 








0.2 







-4.0 



-3.5 -3.0 -2.5 
logio I/roI 



-2.0 




-3.5 -3.0 -2.5 
logio I/roI 



TM (from top): 




^ACDM^ Q 7^ Q Q 5 


''^^ 




= . 


a 


= 0.6 


-4.0 -3.5 -3.0 -2.5 - 


-2.0 


logio I/roI 





FIG. 2: Le]t: Simulation measurements (points) and model predictions (lines) for the peak enhancement of ^hmC?"), i-e., 
(|/_Ro|, o"8^'~"°^), as a function of |/ho1- Note the approximately logarithmic dependence of a A on |/flo|- Middle: The 
peak enhancement of 5hm('") in the abundance-matched case as a function of |/ho| for different values of the power spectrum 
normalization ag'^^^. Right: Same for the threshold- matched case, including the dependence on the scatter a. Effects from 
scatter are negligible in the abundance-matched case. 



typically is the case for models linked to the late-time 
acceleration of the Universe). 

Since the f{R) simulations have only been run for 
one cosmology and a small set of values of |/ho|i we 
use the halo model to interpolate between the simula- 
tion predictions. We have found that the shape g{r) of 
the profile enhancement (see Fig. [T]), when normalized to 
unity at the peak of the enhancement, is independent of 
l/flol to within a few percent for the simulated values of 
l/flol- In the following, we will adopt g(r) measured for 
|/iio| = 10"'^. Hence, we write 

Q(r) = aA(|/«o|,4^°'','^) (18) 

where A is the peak height predicted in the halo model 
as function of |/ho| and a^'^^^ , the cts a ACDM universe 
would have for a given primordial power spectrum ampli- 
tude, and the scatter a. The halo model predictions are 
described in Appendix |^ a is a fudge factor, which is 
determined by matching to Qsimij) at the peak; in other 
words, we are only using the halo model to predict the 
scaling with fno and erg, while the simulations are used 



to match the precise amplitude, a depends on whether 
we are considering the AM or TM case. In the AM case, 
scatter effects on a and g can be neglected, i.e., g is only 
a function of r, and a = 0.52. In the TM case, we have 
a{cT = 0) = 0.73, a{a = 0.6) — 0.77, and interpolate a 
and g{r) linearly in cr. 

For |/flo| and cr^^^^, we use an interpolation based 
on the halo model for |/flo| < 2 x lO^^ and cr^^DM g 
[0.7,0.9]. In order for the MCMC runs to converge, 
however, we need to cover a larger parameter space in 
|/ii:o| and a^'^^^ than can reasonably be covered by 
the halo model. Thus, when |/_ro| > 2 x lO^'^ and 
^ACDM ^ [0.7,0.9], we use an extrapolation fitted to 
the halo model predictions for \fRo\ < 2 x 10~^ and 
^ACDM g jQ Q gj^ described in Appendix[Bl However, 
the details of this extrapolation are not important for the 
final parameter constraints since they lie well within the 
region that is covered by the simulations and the halo 
model inter- and extrapolation (sec WI|) . 

Finally, the prediction for the halo-mass correlation 
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function in f{R) gravity is given by 



(19) 



where here and throughout Chm°^('') is the ACDM 
prediction interpolated from the measurements in the 
ZHORIZON simulations. 



B. Phenomenology with a Gaussian fit 

In addition to the consistent, abundance-matched con- 
straints on /(i?) gravity, we also consider a phenomeno- 
logical approach modeled on the profile enhancement in 
f{R) at fixed halo mass (TM case). This case serves 
to illustrate the ability of halo profiles to probe gravity, 
independently of halo abundances and the specific /(i?) 
model. To do this, we fit Qsim('') for the Ifml = 10~^ 
threshold-matched case without scatter for the ampli- 
tude, width, and position of a Gaussian function in Inr 
and then take the amplitude Fq to be the free parameter 
controlling the modification, i.e.. 



Q™(r,Fo) =FoCxp 



1 / In r — /i 
2 



(20) 



The minimum for the fit of the fixed mass simulation 
(see Fig. [T]) is obtained for e'' = lA7h~^ Mpc and ~ 
l.59h-^ Mpc. Note that Qsim('') in the AM case is not 
simply described by a Gaussian enhancement. 

In the middle panel of Fig. [1] we show the enhance- 
ment of the modified relative to the ACDM (|/i?o| = 0) 
simulated density profile for |/i^o| = 10"'^ and the corre- 
sponding Gaussian function. In the following, we refer to 
this approach as the phenomenological model (PM) case. 

For comparison, Fq matches the peak height of the 
enhancement in the threshold-matched scenario for 

Fo aA{\fRo\ = 10-^, a^^^^ = 0.8,cr O) ~ 0.306. 

(21) 

In general, one can map Fq to the corresponding value 
of I //JO I in the TM case through the right-hand panel 
of Fig. m We shall, however, not restrict the likeli- 
hood analysis to only no n- negative values of Fq, in cor- 
respondence with I //JO I > but extend it to cases where 
Fq < 0, i.e., models where gravity is weakened and pro- 
files are consequently suppressed. A suppression of this 
kind may, for instance, be observed in self-accelerating 
DGP braneworlds [33. 



C. Lensing predictions 



as probes of gravity, in particular, if external information 
on the scatter is available. Note also that the profile en- 
hancement is significantly smaller in the TM case when 
compared to AM at a fixed value of |/flo|- 

We first determine AE^^^m ^^^^-^ i.e., with- 

out including modified gravity effects, using Eqs. (|lip 
through p3p for each concordance model cosmology 
in Table [ill At each r±, we use a four dimensional 
paraboloid to interpolate AS^'^'-"^ in the parameters 
{rim, cr^^'^'^, 71s, cr} . The paraboloid is defined by three 
simulations in each parameter direction. We then inter- 
polate linearly in logr_L. 

In order to include the modified gravity effects for the 
AM/TM case, we write 



AS(ri,|//?o| 



ACDM 



xAS 



ACDM 



(22) 



where AS^'-'^^ contains the dependency on the cosmo- 
logical parameters, Afid 
0.8), and 



^ACDM 



AQfld(r±) 



ASfid(r±,|/flo| ^10' 

AI]^CDM(,^) 



- 1 



(23) 



is obtained by inserting Eq. (|18p into Eq. when per- 
forming the projection, Eqs. pT|) through (|T3)) . using the 
fiducial values for the cosmological parameters defined in 
Table nil with ct = 0. 

Similarly, for the PM case we write 

AS(r^,i^o) - [1 + FoAQ™(r^)] AE^^dm^^^) ^34) 



where 



AQL^(r^) 



AS]fid(ri,Fo = 1) 



AE^CDM, 



(25) 



is obtained by inserting Eq. (|20p into Eq. when per- 
forming the projection, Eqs. (fTTj) through (fT3)) . with fidu- 
cial values for the cosmological parameters. Eqs. and 
()24p are approximate and assume that the r-dependence 
of the modified gravity effects does not depend on the 
cosmological parameters. We found that this approxi- 
mation is valid to better than 1%. 

The effects of varying cosmological parameters on AS 
are illustrated in Fig.|4l Comparing Fig. [Hand Fig.[3l we 
see that the relative enhancement observed in the halo 
profiles in /(i?) gravity is broadened and propagated to 
larger radial scales by the projection and conversion to 
the excess surface mass density. 



Fig. [3] illustrates the effect of varying the cosmological 
parameters and mass scatter a on ^hm- It is apparent 
that the f{R) field strength |/flo| and a have the largest 
impact on the profiles; this shows that halo density pro- 
files in the region of one to a few virial radii are useful 



V. OBSERVATIONS 

The observations in this paper are derived from the 
SDSS [55} , which imaged roughly tt steradians of the 
sky, and followed up approximately one million of the 



8 



^ACDM^ 0.9, 0.7 

n,= 0.95, 1.05 
n,„ = 0.2, 0.3 



— — 



y 



2.0 F 



1.5 - 



— I/roI= 10-^' (AM) 



1.0; 



0.5 



0.0 



0.2 0.5 1.0 2.0 5.0 10.0 20.0 
r [/7 "'Mpc] 



From top: 

o- = 0.4, 0.6, 0.; 



0.2 0.5 1.0 2.0 

r [h "'Mpc] 



5.0 10.0 20.0 







I/roI = 10-^ (TM) 




0- = 




- 0- = 0.6 








1 .. 1 .... 1 1 





2.0 



1.5 - 



- ^ 

1.0 



0.5 - 



0.0 



0.2 



0.5 1.0 



2.0 



5.0 10.0 20.0 



r [/?"'Mpc] 



/ s 
/ \ 
/ \ 



Fo = 1 (PM) 

0.2 0.5 1.0 2.0 

r [h "'Mpc] 



5.0 10.0 20.0 



FIG. 3: Effects on the halo density profile ^hm from varying the cosmological parameters with respect to the fiducial case. 
Upper left: Different parameter values for ag'~'^^ (dashed), Us (dot-dashed), and (dotted). Upper right: \fRo\ ~ 10~^ for 
the abundance-matched case (dashed) and for the fiducial ACDM cosmology with different values of scatter a (dotted). Lower 
left: |/ijo| ~ 10~^ for the threshold-matched case with a = (dashed) and a = 0.6 (dotted) (with the fiducial case corrected 
for scatter). Lower right: Fq = 1 for the phenomenological scenario (see ^TV B[) . 



detected objects spectroscopically j66l - l68| . The imaging 
was carried out by drift-scanning the sky in photometric 
conditions [H, [TOI in five bands (uqriz) (tiI It^ using 
a specially-designed wide-field camera [73]. These imag- 
ing data were used to create the cluster and source cat- 
alogs that we use in this paper. All of the data were 
processed by completely automated pipelines that detect 
and measure photometric properties of objects, and as- 
trometrically calibrate the data [zi^Zil- The SDSS I/II 
imaging surveys were completed with a seventh data re- 
lease |77| , though this work relies as well on an improved 
data reduction pipeline (Photo v5_6) and updated pho- 




that is part of 



tometric calibration (ubercalibration, [Tj 
the eighth data release, from SDSS-III [2 



A. Lens cluster sample 

We use cluster-galaxy lensing measurements around a 
subset of the maxBCG optically detected cluster sam- 
ple from the SDSS, consisting of 5 891 clusters with 
background sources. The parent sample of clusters 
from which our lens sample is derived consists of 13 823 
MaxBCG clusters [38| that arc identified by conccntra- 
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FIG. 4: Effects on the excess surface mass density AS from varying the cosmological parameters with respect to the fiducial 
case. The lensing data has been rebinned for illustrative purposes (cf Fig. [5]) Top left: crl^^DM (hashed), (dot-dashed), 
ftm (dotted). Top right: |/flo| ~ 10"'^ for the abundance- matched case (dashed) and for the fiducial ACDM cosmology with 
different values of scatter a (dotted). Bottom left: Ifm] = 10""^ for the threshold-matched case with a = (dashed) and a — 0.6 
(dotted) (with corresponding values of a used in the fiducial AE for each case). Bottom right: Fo = 1 for the phenomenological 



scenario. 



tions of galaxies in color-position space using 7 500 square 
degrees of imaging data from the SDSS. The entire sam- 
ple is placed into a single redshift slice spanning 0.1 < z < 
0.3 (zeff = 0.23), and a redshift-dependent richness cut 
in iV2oo (the number of red member galaxies above some 
luminosity threshold) is applied to achieve a rcdshift- 
indcpendent number density of fi = 2 x 10~^ft.^^ Mpc'^. 

The maxBCG sample is particularly well suited for our 
study on halo profiles since the BCG is expected to coin- 
cide with the center of its host halo, i.e., the minimum of 
the potential well. If this assumption is perfectly satis- 



fied, then our analysis is simplified since no modeling of 
the mass distribution around satellite galaxies (including 
assumptions about their hosts, cf. (sij ) is required. 

To ensure that this is the case, and reduce effects 
from possible "satellites" (in reality, clumps of galax- 
ies within some larger cluster that are misidentified as 
a separate, nearby cluster) contaminating the maxBCG 
sample, we define a cylindrical region around each clus- 
ter with a transverse radius of three virial radii, derived 
using the mass-richness relation from [82j . and extent 
along the line of sight of Az = ±0.045 (corresponding 
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to X ~ ±100/i~^ Mpc, a 3cr photo-z error). If there is a 
lighter cluster candidate in this region, then the lighter 
cluster is removed from the sample. This removes 30% 
of the clusters in the sample, resulting in a net observed 
number density of n = 1.4 x 10~^h^'^ Mpc'^. As described 
in §111 C[ carrying out the same procedure on the halo 
catalog of the A^-body simulations removes 20% of the 
halos. This finding suggests that of the 30% that were 
removed from the maxBCG sample, 10% were truly spu- 
rious detections and 20% were removed due to chance 
projections. We thus estimate the true parent sample 
number density to be fi = 1.8 x 10~^h~^ Mpc^. This is 
the value used when abundance-matching the halos from 
the f{R) and ACDM simulations. We emphasize that 
it is not a problem that our procedure is overly conser- 
vative; it is better to avoid modeling difficulties at the 
expense of losing 20% of the real clusters in the sample. 
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B. Source catalog 

The catalog of source galaxies (1.18 arcmin^^) with re- 
solved shape measurements and photometric redshifts is 
described in detail by Reyes et al. [s^. In brief, the cor- 
rection for the effects of the point-spread function (PSF) 
uses a method called re-Gaussianization [s^l, with orig- 
inal systematics tests presented in (ssj and an updated 
treatment by Reyes et al. [s^ . The effect of errors in the 
ZEBRA photometric redshifts [s^ on the lensing signal 
calibration was studied by [s^l • 



C. Lensing measurements 

A description of the procedure for calculating the lens- 
ing signal can be found in Reyes et al. [ssj . In brief, we 
assign optimal weights wxs to each lens-source pair based 
on the noise in the shape measurement and based on 
the critical surface mass density Ec''-' = TicntizijZg) esti- 
mated using the source photo- z. To estimate the lensing 
signal AE(r_L), we then compute a weighted average 



AI](rx) 



(Is)^(ls) 



Els wis7r ^ 
27eEr 



(26) 



in logarithmic radial bins. The denominator includes the 
sum over weights of random Zens-source pairs w,s , to cor- 
rect for the dilution of the source sample by "sources" 
that are actually associated with the cluster and arc not 
lensed by it. The factor of 2 arises due to our cUiptic- 
ity definition, and TZ is the shear responsivity, which de- 
scribes how our ellipticity definition responses to a shear 
jssf . After computing the signal, we also compute the 
signal around the random points to check for any sys- 
tematic shear contamination [ssj . and subtract it from 
the real signal (in practice, for the scales of interest, this 
correction is only nonzero for R > 10/i~^Mpc and even 
then, it is well below the statistical errors). Errors are 



FIG. 5: Excess surface mass density, Eq. (|26p . as a func- 
tion of the comoving transverse separation to the cluster cen- 
ter (BCG), r±, measured in the maxBCG sample (points). 
The lines show the predictions from the fiducial and best- 
fit ACDM models as well as for the best-fit phenomenological 
scenario (see t )II B|) . Note that the best-fit abundance- matched 
f{R) model is indistinguishable from the best-fit concordance 
model (see Table [TTT| and is therefore not shown separately. 



calculated using jackknife resampling; for this purpose, 
we divide the survey area and therefore the lens sample 
into 100 equal-area regions. 

We use the same procedures as in Reyes et al. [s^ to 
assess the impact of various sources of calibration biases 
on the lensing signal, and we then remove them, assign- 
ing an overall 5% calibration uncertainty. We therefore 
divide the theoretical predictions for AE by the calibra- 
tion factor C = 1.08, and include a Gaussian scatter of 
0.05 on C when comparing to the lensing measurements 
in the MCMC analysis (see fjlV|. 

Fig. [5] shows the measurement of the unbiased excess 
surface mass density AT,(r±) (multiplied by C) along 
with the best-fit signals for the ACDM, AM, and phe- 
nomenological model (see ^IVB[) . respectively. 

In WI C[ we shall discuss further possible systematics, 
especially those which have scale-dependence. 



D. External priors 

In order to prevent degeneracies of \fRo\ with other 
cosmological parameters and combinations thereof, we 
further employ measurements of the background ex- 
pansion history and the cosmic microwave background. 
For this purpose, we consider the likelihood distribu- 
tion for the concordance model parameters from Q. 
This analysis uses the CMB anisotropy data from 
the five-year Wilkinson Microwave Anisotropy Probe 
(WMAP) [13, the Arcminutc Cosmology Bolometer Ar- 
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ray Receiver fACBAR ) lOOll , the Cosmic Microwave Back- 
ground Imager (CBI) [91|, and the Very Small Sky Array 
(VSA) [ol]. It further utihzes data from the Supernova 
Cosmology Project (SCP) Union [o^l compilation, the 
measurement of the Hubble constant from the Super- 
novae and Ha for the Equation of State (SHOES) [H] 
program generalized by j95| , and the B AO distance mea- 
surements of [9^. For the description of these ob- 
servables, in particular, for the CMB, a high-redshift 
parametrization was chosen, constructed from the physi- 
cal baryon and cold dark matter density fib^^ and Vt^h? , 
the ratio of the sound horizon to angular diameter dis- 
tance at recombination multiplied by 100, 6, the optical 
depth to rcionization r, the scalar tilt rig, and amplitude 
As at K = 0.002 Mpc"^ 



For our analysis we restrict to the parameters that 
are used for predicting the excess surface mass den- 
sity AE in ijlV C[ i.e., Ug and the derived param- 
eters, the total matter density fim and the power 
spectrum normalization cr^^^'^. Hence, we marginal- 
ize over {^]b/l^^^c^^6',T,ln[10l%]} to obtain a three- 
dimensional posterior distribution for rig, f2m, and 
^ACDM ^ -^^jjif.]^ serves as our prior within the MCMC anal- 
ysis. 

Note that by construction, at high redshifts, f{R) 
modifications become negligible, i.e., at large multipoles 
of the CMB, predictions from f{R) gravity match the 
predictions from the concordance model. Modifications 
appear only at low multipoles of the CMB due to the 
Integrated-Sachs Wolfe effect and lead to constraints on 
l/flol of around unity The background expansion his- 
tory within the Hu-Sawicki f{R) gravity model matches 
the one of ACDM for \fm\ <C 1 at the accuracy level 
of current observations. Since we are interested in con- 
straints on f{R) modifications that originate from the 
halo profile alone, we restrict to the concordance model 
predictions for comparison with the data described here. 

As a prior on the scatter a we adopt the probability 
distribution shown in the top panel of Fig. 3 in ob- 
tained from comparing cluster richness with X-ray mass 
measurements. This constrains the scatter to be < 0.7 at 
the 95% confidence level. While that analysis assumed 
GR, the measurement of the scatter in the mass-richness 
relation only relies on the fact that the X-ray mass prox- 
ies trace true mass with much smaller scatter than rich- 
ness. This is expected to hold even in the modified grav- 
ity case, at least when the chameleon mechanism is not 
active [93| as is the case for the values of | //^o | considered 
here. 



Finally, for the lensing calibration, which we use to 
scale AE (see S|V|, we use a Gaussian distribution around 
1.08 with 5% standard deviation. 



VI. RESULTS 

We now move to the MCMC likelihood analysis of the 
cosmological parameter spaces 



and, in the case of the PM enhancement. 



ACDM 



3, cr,C, Fq} , 



(27) 



(28) 



where for the concordance model Pacdm = Vam H 
= 0} = VpM n {Fo =0}. We implement the fol- 
lowing flat priors on the parameters in 7^am\'Pacdm and 
7'pm\7'acdm: \fRo\ G (0,10) and F^ € (-5,5) for the 
AM and PM enhancement, respectively. In addition to 
the priors from the distance and CMB measurements dis- 
cussed in WD[ we further employ flat priors on top of 
the priors on the parameters in T^acdm^ ^^m S (0.05, 0.5), 
^ACDM g (0.4,1.6), ris G (0.5,1.5), a € (0,2), and 
C £ (0.5, 1.5). Note that these bounds only serve as clear 
truncations for the parameter exploration in the MCMC 
code and since the ranges are chosen much wider than the 
bounds from the external priors in ijVDI and of C in 
they do not affect the final parameter constraints. 

The COSMOMC [ill package used for the MCMC like- 
lihood ana lysis employs the Metropolis-Hastings algo- 
rithm [9^, llOOf fo r th e sampling and the Gelman and 
Rubin statistic G |lOl| for testing the convergence. We 
require — 1 < 7 x 10""^ for our runs. We summarize our 
results in Table iHll 



A. f{R) gravity 

Fig. |4] shows that the f{R) predictions in the AM case 
for I /ho I = 10~'^ are in clear tension with the data, 
at least when excluding scatter. Further, due to the 
strong radial dependence, the f{R) effects cannot easily 
be canceled by varying some of the other cosmological 
parameters. This leads to a ID-marginalized constraint 
of Ifml < 3.5 X 10~^ at the 95% confidence level. Note 
that including a prior on scatter plays an essential role 
(see Fig. [B]), i.e., if we were to remove it from the anal- 
ysis, very large scatter would make large |/i?o| models 
viable (see Fig. [3] and Fig. 2]) and due to a slow increase 
of the enhancement A as function of \fRo\ (see Fig. [2|), 
there would be a rather loose constraint on Ifml- This 
is what happens if one wishes to constrain the TM sce- 
nario instead. In contrast to the AM case, when fixing 
the mass range equally in the ACDM and modified grav- 
ity model and therefore lowering the average halo mass 
within the stacked profiles, the discrepancy in the en- 
hancement on AS on scales below r± < lh~^ Mpc and 
above r± > 10h~^ Mpc is less severe (see Fig. [T] and 
Fig. 21). Therefore, in that case, smaller values of scatter 
are already sufficient to make the profile enhancement 
compatible with the cluster-galaxy lensing data. This 
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Parameter 


ACDM 


AM 


PM 


ACDM 

a 

lO-'l/flol 
Fo 
C 


0.266 ±0.011 0.268 
0.795 ±0.016 0.791 
0.956 ±0.011 0.951 
0.46 ±0.10 0.46 

1.083 ±0.048 1.089 


0.251 ±0.013 0.265 
0.769 ± 0.022 0.788 
0.961 ±0.015 0.952 
0.53 ±0.13 0.45 
< 3.55 0.00 

1.114 ±0.052 1.085 


0.261 ±0.011 0.258 
0.785 ±0.017 0.776 
0.956 ±0.012 0.952 
0.47 ±0.10 0.42 

0.34 ± 0.20 0.34 
1.092 ±0.049 1.084 


-21nL 


14.2 


14.2 


11.5 



TABLE III: Mean, standard deviations, and best-fit values for tlie concordance model, f{R) gravity in the abundance-matched 
case, and the phenomenological model, respectively. For |/flo| we quote 95% ID-marginalized confidence levels. — 21nL is 



calculated for the cluster-galaxy lensing data including the priors of ^VD 
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FIG. 6: 2D-marginalized contour plots for the abundance-matched (top row) and the phenomenological enhancement case 
(bottom row), showing 68%, 95%, and 99% confidence levels. The dashed line corresponds to ACDM predictions from DMO 
simulations. 



leaves |/flo| unconstrained within the region of applica- 
bility of our linearized /(i?) equations, l/^ol ^ 2 x 10~^. 
Note, however, that the TM case does not consistently 
take into account the enhanced abundance of clusters of 
f{R) gravity and we, therefore, restrict our f{R) specific 
constraints to the AM case. 

Fig. |6] shows the 2D-marginalizcd likelihoods for the 
parameter degeneracies with |/i?o| and Fig. [7] shows 
the one-tail ID- marginalized likelihood for |/i?o|- In 
Fig. m we illustrate the band of ^hm and AS predictions 
bounded from below by the best- fit AM f{R) gravity 
model, which is essentially identical to the best-fit con- 
cordance model (l/flol = 1-7 x 10~^), and from above 



by the upper 68% confidence level value of |/flo| with 
otherwise identical parameter values. 



B. Phenomenological scenario 

In the phenomenological case based on the TM halo 
profile enhancements, we use a Gaussian function in Inr 
with width and position fixed to fit the simulation and 
only consider the amplitude of the Gaussian function _Fo 
as an additional free parameter (see ijIVB[) . We obtain a 
mean and standard deviation of Fq = 0.34 ± 0.20. Fig. [S] 
illustrates parameter correlations with Fq, Fig. [7] shows 
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I/roI ^"0 



FIG. 7: One- and two-tail ID-marginalized likelihood. The dotted lines indicate the 68%, 95%, and 99% confidence levels, the 
dashed line corresponds to the ACDM prediction from DM0 simulations. Left: \fRo\ in the abundance-matched case. Right: 
Fo in the phenomenological scenario with a Gaussian fit in Inr to the enhancement in the TM case. 
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FIG. 8: Left: Best-fit prediction for ^hm with respect to the best-fit ACDM model prediction for the phenomenological scenario 
(dashed) and fiducial ACDM cosmology (dotted), respectively. Right: Best-fit prediction for the excess surface mass density 
AE in the phenomenological scenario (dashed) and the fiducial ACDM model (dotted) with respect to the best-fit concordance 
model. Note that the lensing data has been rebinned for illustrative purposes (cf. Fig. [5]). The shaded areas indicate regions in 
the abundance-matched case bounded by the best-fit model from below, and |/ho| = 3.5 x 10"'^ (corresponding to the 68% CL 
bound) with otherwise identical parameter values from above. The best-fit model for the abundance-matched case is essentially 
identical to its ACDM counterpart and is therefore not shown separately. 



the two-tail ID-marginalized likelihood for the amplitude 
Fq, and in Fig. [5]wc present the best-fit predictions for 
^hm and AS. The best-fit parameter values, as well as 
the corresponding — 21nL are listed in Table Hill 

Fq = corresponds to the ACDM model with DM0 
simulations and is consistent at the ID-marginalizcd 95% 
confidence level (Fig. O. The best fit is obtained for 



Fo = 0.34, achieving — 2AlnL = —2.7 with respect to 
the best-fit concordance model. 

The data thus slightly prefer an enhancement in halo 
profiles over ACDM in the phenomenological case. Fu- 
ture surveys will thus cither strengthen the constraints 
on modified gravity parameters, or even more interest- 
ingly, provide additional evidence for Fq > 0. Note that 
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for the best-fit ACDM model the reduced is roughly 
unity and that we therefore do not expect our error bars 
to be significantly underestimated. 

C. Systematic effects 

The shape of the enhancement effect, g{r), on the clus- 
ter profile ^hm and the excess surface mass density AS 
observed in f{R) gravity simulations cannot be repro- 
duced by any reasonable deviations in the parameter val- 
ues of the fiducial cosmology (see Figs. [3]and|4|). Our 
comparison of theoretical predictions to the lensing ob- 
servable is, however, affected by the following possible 
scale-dependent systematics. 

• Mass scatter: In mil C[ we include a log-normal 
scatter in the mass-richness relation in our theo- 
retical modeling. The actual form of the scatter 
might differ from log-normal, though it seems un- 
likely that this would result in an enhancement of 
AE localized at ~ (1 - 10)11^^ Mpc. 

• Baryons: In order to understand the formation of 
galaxies within clusters, it is essential to include the 
baryonic components. Realistic models comprise 
mechanisms such as gas cooling, star formation, su- 
pernovae feedback, as well as the feedback from su- 
permassive black holes to avoid the overcooling and 
accumulation of gas in the core of the cluster, the 
so-called active galactic nucleus (AGN) feedback. 
AGN outbursts produce shock waves that move the 
gas from the core to larg er ra dii, i.e., between 
and 2rv, as was shown in jl02| by employing simu- 
lations of Virgo-like galaxy clusters. Moreover, due 
to the AGN feedback, there is a slight adiabatic 
expansion of the dark matter when compared to 
DM0 simulations (see Fig. 9 of |l02| ) , leading to a 
< 10% effect on the density profiles. These effects 
have a radial dependence that is qualitatively differ- 
ent from the modified gravity enhancements consid- 
ered here. Note that the moved mass by modified 
gravity can be calculated as 

/•CJO 

AM = Attp ^u„,{r)Q{rydr (29) 
Jo 

and amounts to AM w 6 x lO^^h^^ Mq for the 
phenomenological fit with _Fo ~ 0.3. 

• Intrinsic alignment: High-precision weak- lensing 
measurements may be contaminated by th e intrin- 
sic alignment of galaxies (see, e.g., |103| ). The 
correlation of intrinsic alignment and gravitational 
shear distortion can contribute to the observed el- 
liptic ity correla tion function and AS at the < 10% 
level pMflOSj . 

• Miscentering and satellites: The cluster centers in 
the MaxBCG sample are identified by the bright- 
est cluster galaxy (BCG). The true cluster center 



may, however, be offset from the BCG position (see, 
e.g., discussion in [13). This effect causes a sup- 
pression of the lensing signal in the inner parts of 
the halo, which subsequently leads to an underesti- 
mation of the cluster mass and the concentration. 
A miscentered AS can have a bump relative to a 
correctly centered AS, which is, however, located 
furt her i nwards than the f{R) gravity enhancement 
(cf. |l06l |). A similar enhancement around the virial 
radius can further be introduced by galaxy satel- 
lites. To prevent the contamination of the excess 
surface mass density through satellites, we apply a 
cylindrical cut in the projected radius at rcut = 
in the simulations (see ^III Bp and the observations 
(see i fV|) . Note that we applied this cut only to the 
ZHORIZON simulations and not to the f{R) grav- 
ity simulations. We verified, however, that chang- 
ing Tcut has a negligible impact on the dependence 
of AS on cosmological parameters. Furthermore, 
the cut only affects AS on scales r± > 5h^^ Mpc. 
Therefore, we can safely assume that there is no sig- 
nificant impact on the relative f{R) enhancement 
by the cylindrical cut. 

• Wrong cosmology: The analysis of lensing as used 
in this study requires the assumption of an a pri- 
ori cosmological model to estimate the critical sur- 
face mass density Scrit and to convert angles to dis- 
tances. Within ACDM, a wrong prior on the cos- 
mological model produces a radial horizontal shift 
of AS at the < 2% level for rj,„ = 0.25 ± 0.05 (see 
discussion in jl07( ). Note that the Hu-Sawicki f {R) 
gravity model matches the ACDM background to 
order |/i?o|- Deviations of this magnitude have a 
negligible impact on AS. 

• Simulation systematics: In order to test the conver- 
gence of the halo profiles of the large-scale cosmo- 
logical simulations on the scales used in this study, 
we compared the halo profiles from the zhorizon 
simulations to th e halo profiles of the millennium 
simulations jlOSj . which employ N = 2 160'^ parti- 
cles in a 500'^h~^ Mpc^ box. The profiles agree at 
the 5% level on the scales of interest. We there- 
fore conclude that the zhorizon simulations have 
converged for r ^ (0.2 — 100)h~^ Mpc. The halos 
in the zhorizon simulations are identified using an 
FoF halo finder, while the f{R) effects were mea- 
sured on a SO-identified halo sample. However, 
since we only use the enhancements from the f{R) 
simulations relative to |/i?o| = 0, we expect the dif- 
ference to be smaller than the residual statistical er- 
ror 20%) on the modified gravity effects. Note 
that the environmental effects found in [93, Il09j 
are induced by the chameleon mechanism and are 
not relevant for the values of |/ijo| considered here. 
Moreover, halo finders typically agree at t he sc ales 
relevant to our halo profile measurements jllOl |. 

• Survey geometry: While we are mimicking the se- 
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FIG. 9: Current constraints on f{R) gravity. On linear scales, 
the strongest bound on \fRo \ is obtained through the compar- 
ison of predicted to observed cross correlations of the ISW 
with foreground galaxies. In the nonlinear regime, enhance- 
ments of the abundance of clusters and the cluster density 
profile due to the f{R) modification are incompatible with 
observations unless |/i?o| is smaller than 10~* and 10"'^, re- 
spectively. The currently strongest bounds on |/flo|, however, 
are inferred from requiring the modification to be suppressed 
by the chameleon mechanism within the solar system and the 
dark matter halo as well as from strong gravitational lenses. 



lection process as closely as possible, including the 
removal of fake clusters, the simulation measure- 
ments provide the dark matter and halo positions in 
a cubic box. Furthermore the simulation results are 
obtained from a single redshift slice at the effective 
redshift of the sample, and are thus not accounting 
for the redshift evolution of the lens sample. 



Except for the case of the scatter, we neither model 
the systematics described above nor include them as ad- 
ditional errors to the measurement when performing the 
likelihood analysis. In order to consistently include these 
systematics, they should not only be carefully analyzed 
within ACDM but also in the context of f{R) gravity, 
which is beyond the scope of this paper. Note that, 
when added in quadrature, the described uncertainties 
sum up to a < 15% and < 25% error in the predicted 
AS for ACDM and the modified gravity cases, respec- 
tively. This work is based on the assumption that the 
above systematics, except for the mass scatter, can be 
neglected and that the observations can correctly be de- 
scribed by an average over the DM0 simulations. Note 
that our ACDM model indeed provides a good fit to the 
data (see Fig. [S]). 



VII. CONCLUSION 

Modifications of GR as in the f(R) gravity model un- 
der consideration in this paper generically predict depar- 
tures from the standard growth produced in the concor- 
dance model. On the largest, cosmological scales (r > 
10 Mpc) and on small, solar-system scales (r < 20 AU) 
such deviations have extensively been instrumentalized 
to probe gravity. However, structures on intermediate 
scales also offer opportunities to test the gravitational 
interactions. 

In this paper, we test modifications of gravity on 
scales around the virial radius of a cluster, i.e., r ~ 
(0.2 — 20) Mpc. The modification of the Poisson equa- 
tion leads to a difference in the accretion of mass onto 
massive dark matter halos. The resulting halos exhibit 
enhanced density profiles at a few virial radii that of- 
fer a unique opportunity for testing gravity. We use the 
projected mass distribution measured through cluster- 
galaxy lensing around maxBCG clusters from the SDSS 
to put constraints on the modifications induced by the 
Hu-Sawicki f{R) gravity model. For consistent theoreti- 
cal predictions we rely on f{R) gravity and concordance 
model A^-body DM0 simulations. Matching simulated 
to observed halos by abundance, we obtain a one-tail up- 
per bound of l/flol < 3.5 x 10"^ at the ID-marginalized 
95% confidence level. This places a new independent 
constraint on f{R) gravity at intermediate scales, where 
\fm\ < few 10-4 and |/flo| < (10^^ - lO'^) are current 
bounds inferred from large cosmological and solar-system 
scales, respectively. We summarize current constraints 
on Ifml in Fig. 1^1 showing bounds inferred from measure- 
ments in the solar system JH, of strong lenses [IJ], the 
abundance of clusters 0, la7galaxy-ISW (gISW) cross 
correlations 0, 0| , and the CMB 0, 0] , as well as our 
constraint from halo density profiles measured via weak 
gravitational lensing. We extrapolate results presented 
in to estimate an upper bound on |/jjo| from the Eq 
measurement of Q, which combines weak lensing mea- 
surements around galaxies with their large-scale veloci- 
ties. Note that Fig. [3] does not include the measurement 
of gravitational redshifts of galaxies in clusters at around 
(1 — 6)h~^ Mpc of [l^ since it was found to be consistent 
with f{R) gravity and cannot be illustrated in the same 
manner as the previous measurements. 

In order to assess the ability of halo profiles to pro- 
vide constraints independently of halo abundances, we 
also considered a phcnomenological parametrization of 
the modified gravity effects on halo profiles at fixed mass. 
In this scenario, the concordance model (with amplitude 
of the modification Fq = 0) is consistent with the lens- 
ing measurement at the 95% ID-marginalized confidence 
level. Thereby, we considered a Gaussian enhancement of 
the cluster density profile due to modified gravity located 
at a few virial radii. The best-fit value of Fq = 0.34±0.20 
indicates that the data slightly prefer an enhancement in 
halo profiles over ACDM; future surveys will thus either 
strengthen the constraints on modified gravity parame- 
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ters, or even more interestingly, provide additional evi- 
dence for Fq > 0. 
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Appendix A: Halo model predictions for the density 
profiles 

In this appendix, we describe the halo model prediction 
for Q(r), the enhancement in ^hm('^) induced by f{R) 
gravity, Eq. (|18p . In the halo model, the halo-mass cross- 
correlation function is given by a sum of two terms. 



2h 



(Al) 



denoting the 1-halo and 2-halo contributions, respec- 
tively. Throughout, all quantities are evaluated at the 
redshift of the f{R) simulation output, z — 0.22. Por the 
TM case, we consider halos with mass Ma > Mq, where 
A = 300 is the overdensity in units of the background 
matter density today, and Mq = 10^^-^^h~^ Mq is a fixed 
threshold mass determined by matching the halo abun- 
dance in simulations to the observed abundance, i.e., the 
same for f{R) and GR. In the AM case, Mq is determined 
separately for f{R) and GR through 



rivdlnil/v 



(A2) 



In J\/v,o 



where is the mass function of dark matter halos per 
logarithmic interval in the virial mass My = Ma„. Wc 
adopt a fixed virial overdensity of Av = 390. The virial 
mass threshold is then converted to the threshol d M q 
for A = 300 through the rescaling described in |lllj . 
For the viria l mas s function, we adopt the Sheth- Tormen 
prescription [ll2j |. 



dn 



d In Mv 



where = 6c/cr{M), a{M) being the variance of the 
density field for a top-hat enclosing mass M at the back- 
ground density, and 

vf{iy) = A\j'^av-^ [1 + {au^)-P] e-'^'^'/2 (A4) 

with a = 0.75, p = 0.3, and Sc = 1.673. A is fixed so that 
Sdvf{v)^l. 

The 2-halo contribution to ^hm('') is most easily writ- 
ten in terms of its Fourier-space counterpart P^m{k), de- 
fined through 



" p2h(i \ ik-r 



(A5) 



The two-halo halo-mass power spectrum is given by 

pZ = bi> Mn)I{k)P,n{k), (A6) 
where P,n is the linear matter power spectrum, 
lZj.J{My)nAM.)d\nMy 



b{> Mo) 

m 



S^AU,,MM.)d\nMy 
'° M 

n^^yik, My)b{My)dln M^ , 

Pni 



(A7) 
(A8) 



and b{Mv) is the scale-independent linear peak- 
background split bias derived from the Sheth- Tormen 
mass function: 



6(A'/v) = 6(A: = 0,Mv) 



Sc 



(5c [1 + (az/2)P] 



(A9) 



y(k, M) is the Fourier transform of a Navarro-Prcnk- 
White (NFW) [m density profile which is truncated 
at the virial radius r^. y is normalized so that y{k = 
0, M ) = 1. We adopt the mass-concentration relation of 

The one-halo contribution is simply the normalized 
stacked NFW profile, 

?hmW = A^"' / /ONFw(r; A/v) d\nM, (AlO) 

Jin AU 

N= I nvdlnAfv. (All) 

Jin A/v 

The halo model prediction for ^hm(^) in f{R) gravity 
is then obtained by substituting the linear f{R) matter 
power spectrum into the above expressions. Wc do not 
change the concentration relation, motivated by the fact 
that the inner profiles o f halos seem relatively little af- 
fected by f[R) Ell- The resuh is sho wn in Pig. [TOl 
scaled by the factor a introduced in i jlll C[ together with 
the simulation results. The halo model prediction pro- 
duces a bump at a few virial radii, because halos arc on 
average more massive in f{R) gravity, leading to slightly 
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FIG. 10: Left: Halo model prediction scaled by the overall factor a — 0.52 (Sec. IIII C|) in comparison with simulation 
measurements for the abundance-matched case and \fRo \ = 10"^ Middle: Same as left panel, but for the threshold-matched 
case without scatter (halo model scaled by a = 0.73). Right: Same as middle panel, but with a scatter of cr = 0.6 (halo model 
scaled by a = 0.77). 



larger virial radii. This becomes noticeable because the 
stacked truncated profile becomes very steep outside the 
virial radius corresponding to Mv,o- Clearly, there are 
discrepancies between the halo model predictions and the 
simulation results at both small and large r. Hence, we 
only rely on the halo model prediction for the overall 
amplitude, whose scaling as function of |/i^o| is predicted 
well (Fig.m, whereas the radial dependence is taken from 
the simulation measurements. 

As an aside, in (36| . modified spherical collapse param- 
eters were derived for a collapse with enhanced forces 
throughout (i.e., the limiting case of infinite reach of the 
fifth force). This set of parameters can be used to esti- 
mate the spread in the halo model predictions induced 
by modified gravitational forces. We found that the halo 
profile predictions for both sets of spherical collapse pa- 
rameters are very similar, with the unmodified param- 
eters yielding a somewhat better approximation to the 
simulation results. 



Appendix B: Extrapolation of the halo model 
predictions 

In order for the MCMC runs to converge, the chains 
need to cover a large parameter space. At the most ex- 



treme values of the cosmological parameters, the halo 
model approach (see i ^IV A[) breaks down and we need to 
rely on a more ad hoc extrapolation for A ( | /j^g \,o's'^^^)- 
We design it to fit the simulations and halo model pre- 
dictions within the AM scenario in the range l/^ol < 
2 X 10^2 and cr|^c:DM ^ jq g gj ^gjj^g ^j^g functional form 

A ^ ao{as) + ai{as) X + a2icrs) e'^ , (Bl) 

where x = log^o |/iio|- The approximation, Eq. (|B1[) . is 
accurate at the < 0.1% level within the range of simu- 
lated values of Ifml, i.e., the regime of correspondence 
to the f{R) gravity model. For the coefficients, we use 
the fit 

ai{as) = ttio + anas + ai2al, (B2) 

where i = 0,1,2. Note that we used as = a^'^^^ in 
Eqs. (|B1[) and (jB2[) to simplify notation. 

Furthermore, note that the exact form of this extrap- 
olation does not affect our constraints, which are well 
within the halo model inter-/extrapolation. 
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